Mean-field theory for Bose-Hubbard Model under a magnetic field 
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We consider the superfluid-insulator transition for cold bosons under an effective magnetic field. 
We investigate how the applied magnetic field affects the Mott transition within mean-field theory 
5^ ■ and find that the critical hopping strength (t/U) c , increases with the applied field. The increase 

in the critical hopping follows the bandwidth of the Hofstadter butterfly at the given value of the 
magnetic field. We also calculate the magnetization and superfiuid density within mean-field theory. 

Q\ • PACS numbers: 05.30.Jp, 05.70.Fh, 67.40.Db 
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I. INTRODUCTION 



One of the most interesting developments in ultra-cold atom physics is the study of neutral atoms in optical lattices 
An optical lattice is prepared by creating a periodic potential utilizing a standing wave of light, and optical 
lattices in one, two and three dimensions have been realized experimentally. 

The cooling, trapping and coherent manipulation of the atomic motion by their interaction with light has been 
established by numerous investigations in the field of atom interferometers [|J , matter- wave superradiance Q , matter 
wave parametric amplifiers 0, Q and others 0] . One should mention also the opportunities garnered by using the 
ultra-cold alkali atoms as quantum computers Q, or by using the Mott insulating state of neutral bosonic atoms 
for detection of quantum entanglement [g, IToj .The custom-made trapping potentials in the optical lattice has also 
opened a venue to study many condensed matter problems, such as the Mott insulator transition experimentally 
realized by Greiner et al. [HI 

Although many different regimes exist for optical lattice experiments, one that is quite interesting from a theoretical 
point of view is that of a deep lattice with few particles per lattice site. If fermions are used instead of bosons, these 
experiments may lead to direct realization of many correlated electron model Hamiltonians such as the Hubbard 
model of high temperature superconductivity or lattice Quantum Hall models (T3|. 

In this work, we concentrate on bosons, and assume that at each lattice site there is only one available state (that 
' is equivalent to requiring the first excited state at each lattice site to be sufficiently high in energy). In this case the 
qq , Hamiltonian is(l3j 

o, r _ 

q| H = -t^(a\a j + a\a i ) + —^(n i -l)rH-li^rH (1) 



where ai is the annihilation operator at site i and rii = a\ai is the number operator at site i. The first term corresponds 
to hopping between different lattice sites and for practical purposes only nearest neighbor hopping is important, so 
the sum < ij > is carried over the nearest neighbors. Second term is the particle-particle interaction and the last 
term is the chemical potential. This is the widely studied Bose-Hubbard Hamiltonian [3, [IJ, [l5| . 

The strong tunneling limit between optical lattice sites (U/t << 1) corresponds to the superfiuid (SF) phase. 
Changing the laser intensity with increasing depth of the optical potential the atomic waves become more localized 
and the on-site interaction U increases at the same time with reduction of the tunneling parameter t [TH ] . The system 
is driven to a Mott insulator (MI) phase and loses long range phase coherence. In general, if the interaction is strong 
enough the system prefers a particle number that is commensurate with the number of lattice sites and the system 
goes into the insulating phase. The strength of interaction U needed for this transition is roughly the bandwidth of 
the noninteracting system 2zt, where z is the number of nearest neighbors. 

A much less studied problem is that of the Bose-Hubbard Hamiltonian under a magnetic field. Experimentally, 
of course, the bosons used in cold gas experiments are uncharged and would not be directly affected by an external 
magnetic field. However, recent studies have shown in detail how a magnetic Hamiltonian , or in general effective 
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electromagnetic fields can be generated for atoms in optical lattices using an external time varying electric field flq]. 
by an oscillating quadrupole potential together with a periodic modulation of the tunnelling between lattice sites [171 ] 
or using more complicated laser configurations (H)]. These investigations suggest that it may be possible to study 
such systems with reasonable improvements on already functioning experiments. An effective magnetic field can also 
be created by rotating the optical lattice, and cancelling the centrifugal force of the rotation by an external quadratic 
trap We also note the two recent papers which study how the artificial external non-Abelian gauge potentials 
can be created for cold atom systems [20, [2l| . 

In this work, we assume that we have a two dimensional square lattice in the x — y plane, under a magnetic field 
in the z direction. We also consider that the "charged bosons" are interacting only when they are on the same lattice 
site and the temperature of the system is set to zero. In this case the Hamiltonian is 

W = -t [«L tt Ki)m + e l27ripn al m a n{m+1) + h.c] (2) 

n,m 

+ "^" ^ ^ { a nm a nrn ~~ ^) a nm a nrn ~ A* ^ , a nm a "nm ■ (3) 

Here, we label every site of the lattice i, j by two integers i = (nj,mi) = (n,m), one integer (n) along the x axis, 
the other (m) along the y axis, and choose the gauge to be A = xy. The first term is the usual hopping term, where 
hopping along the y axis gets a phase shift due to the presence of the magnetic field. Magnetic field affects the system 
through the parameter ip with 

V = Bf/ipo (4) 

where / is the lattice spacing and ipa is the flux quantum. Thus, the parameter <p measures the magnetic flux per 
unit cell of the lattice in units of flux quantum . Second and third terms are interactions and chemical potential, 
respectively. Even at U = 0, the noninteracting limit of this Hamiltonian shows interesting results; the energy 
spectrum at U = is known as the Hofstadter butterfly [23] • Most important aspect of the non- interacting problem 
is that the bandwidth depends critically on ip, and gaps open up or close in a self-similar manner. With such a 
complicated single particle spectrum, it is not at all clear how the presence of the magnetic field will change the Mott 
transition, or whether more exotic phases can be found. 

We believe that with the possibility of experimental realization, it is of importance to study this model more 
closely and understand its rich phase diagram. In this work, we concentrate on the superfluid-insulator transition and 
investigate the effect of the external magnetic field on the phase boundary. Our mean-field approach is not capable 
of capturing possible correlated phases, however it may serve as a basis for more detailed investigation of the model. 

We find that the Mott insulating phases become more stable under the applied magnetic field, an expected effect 
as one of the most important effects of the magnetic field would be to localize particles further. More importantly, 
we find that the critical hopping to interaction ratio t/U, roughly follows the bandwidth of the Hofstadter butterfly 
We also calculate the magnetization and the superfluid density within mean-field theory. 

In the rest of this paper we first outline our calculational scheme of solving the Bose-Hubbard model under a 
magnetic field within the mean-field approach. We then present our results on the phase diagram identifying the 
superfluid and insulating regions. We conclude with a brief summary of our main results. 



II. MEAN-FIELD APPROACH 



Our calculations are based on the mean-field approach of the Bose-Hubbard Hamiltonian [23j, by considering the 
following decoupling formula for the product of the two Bose field operators: 



a nm a (n+l)m ~ ( a nm) a (n+l)m + a nm( a (n+l)m) ~ ( a nm) ( a (n+l)m) ■ (5) 

The average value (o4 m ) represents the order parameter ^ nm that accounts for the insulator-superfluid transition. 
It is equal to zero on the insulator side of the transition when the ground state of the system has a definite particle 
number on every site of the lattice, and has a nonzero value for the superfluid state when there are large quantum 
fluctuations of the atom number in the optical lattice. In this case |\& nm | 2 represents the local density of the atoms 
in the condensate state. 
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Using Eq. ([5]) the Bose-Hubbard Hamiltonian given in Eq. (|2|) turns into a sum of the following single-site terms: 

~\~U {jl nra l)71nm [ATljim -\- C nra 3 (6) 



*(li+l)m"nt71 ' T {n-l)m u nm ' c ^ n(m+l) "nm ' L *n(m-l)"ran ' 1J -"-- 



where n nm is the single-site density operator tt£ m ffinm and C nm is a constant energy term. 

The matrix elements of the mean- field Hamiltonian Ti^m m the occupation number base \N nm ) are given by: 

(N nm \7inm \N nm ) — -^UN nm (N nm — 1) — pN nm + C nm , (7) 



(N nm + 1 \H^\N nm ) = - ty/N nm + 1 ** n+1)m + ^ {n _ l)m + e^ n K {m+ i) + e- i2nvn K(m-l) > (8) 



where we used the property of the Bose field operators c\N + 1) = V^V + 1| A?") and c^|iV) = V^V + l|iV + 1). All other 
matrix elements are zero, except the conjugate elements of Eq. ([5]). We note that the occupation number N nm above 
varies from to oo and they are referred to the site (nm) of the optical lattice. We diagonalize the Hamiltonian Eq.© 
in a truncated basis \N nm ) with N nm = O...N max and calculate the ground state of the mean-field Hamiltonian: 

\G nm ) = £ aTWnm) , (9) 
N=0 

with the coefficients a^ m corresponding to the lowest eigenvalue of the matrix of Eq. (|6|) in the truncated base. 
The order parameter corresponding to the ground state given by Eq. © will be 

JV mQX -l 

*nm = (G nm \aUG nm ) = £ aT*aT +1 VN+T. (10) 

N=0 

For a given truncated basis the equations of the finite Hermitian matrix in Eq. ([7]) and Eq. ||5J) and the formula for 
the SF order parameter Eq. ([TO)) represent a set of self-consistent equations that give the solution of the ground state 
of the single site Hamiltonian Eq. ^ and the order parameters ^ nm in the mean-field approximation. 

The numerical calculations are repeated with increasing values of the dimension N max of the truncated basis to attain 
convergence of the solution. In the mean-field approximation, the ground state of the Bose-Hubbard Hamiltonian 
Eq. is given by the direct product of the single site ground states of Eq. 



|G) = -Q|G" m ). 

nm 

For a given ground state Eq. ^ the probability for the single site operator n nm to take the value N will be given 
by the square of the corresponding developing coefficient a™ 1 . The average single site occupation number denoted 
with p(nm) is equal to: 

p(nm) = (n nm )= £ \a n N m \ 2 N, (12) 

and the condensate component of the superfluid density, within mean-field theory, on the site nm is: 

p s {nm) = |* nm | 2 . (13) 

We denote by p and p s the surface average of p(nm) and p s (nm) respectively. For a noninteracting system at zero 
temperature all of the bosons are in the condensate, i.e. in the lowest single particle state of the lattice and we have 
p = p s . When the interaction increases (nonzero values of U in Eq. ^) only a significant fraction of the bosons will 
condense in the same single particle quantum state, and we have p s < p. The competition between the kinetic energy 
of the system t and interaction U gives rise to interesting successive transitions between a superfluid and a Mott 
insulator. 

The mean- field solution of the nonmagnetic system and the phase diagram are calculated by Sheshadri et al. [23[ . 
See also the Mott insulator lobes in Fig. [1] and [5] For the magnetic Hamiltonian, perturbative techniques are used by 
Niemeyer et al. [24[, where the MI lobes are calculated for small values of the magnetic flux <p = 0, 0.125. 
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Introducing the magnetic field in the hopping term of the Hamiltonian of Eq. breaks the temporal invariance 
of the Hamiltonian and gives rise to persistent current flow of the " charged bosons" . For any bond connecting 
neighboring sites (i;k) = (n,m;n± l,m) (or (i;k) — (n,m;n,m ± 1)) of the lattice, we define tu- — t for hopping 
along x (or t^k — te %2lxnip for hopping along y) . We calculate the local bond current of the superfluid phase using the 
following formula: 



1 r 



in 



tik<4 a k ~ tua\ai , (14) 



by substituting (a|afc) with (^VE^). Another parameter of interest is the magnetic momentum. In the mean-field 
decoupling we can define the single site magnetization by the following formula: 

Mnm — nv x (nm) — mv y (nm) = — ^p- (15) 

Trn dcp 

where fi nm is the average value of the mean-field Hamiltonian of Eq. ([6]) with respect to the ground state Eq. ([9]) and 
represents the single site energy of the Bose gas. The averaged site velocities v x ' y (nm) are equal to the average of the 
bond currents of Eq. (j 14(1 connecting the site n, m to its neighbors along x or y accordingly. The magnetic momentum 
denoted with M is equal to the surface average of Eq. (TT5")) . 

For zero magnetic field the 2D lattice has the translational invariance along both axes and the order parameter 
$n m does not depend on the site index (nm). 

In our case, for nonzero magnetic field in the chosen Landau gauge the system preserves only the invariance along the 
y axis of the lattice. Therefore, the order parameter is chosen as ^ nm = ^>n- I n this case, we calculate the mean-field 
solution for a given ratio of the magnetic flux <f> = p/q. From the equation of the matrix elements of the mean-field 
Eq. ©, it can be noted that the periodicity of the solution is ^„ = ^ n + q . The same periodicity condition is verified 
by the the density p that also shows the translational invariance along the y axis: p(nm) = p(n) and p(n) = p(n + q). 
To solve the mean-field equations we choose a finite sequence of the lattice of dimension q in ^-direction and impose 
periodic boundary conditions. In the Landau gauge, the system is periodic with lattice periodicity in the y direction, 
thus our calculations are carried out on a 1 x q supercell with periodic boundary conditions. We measure the energy 
in units of U (i.e. U = 1). 



III. RESULTS AND DISCUSSION 



Based on our calculations outlined in the preceding section, we now present our numerical results. The two possible 
states of the 2D Bose system are selected as follows: in the MI phase, the on-site occupation number has integer 
values and the variance of p is zero; in the SF phase the on site occupation number has noninteger values and the 
variance cr(p) ^ 0. The main features of our results are illustrated in Figs. 1-6. 
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FIG. 2: (color online) Phase diagram for the magnetic flux ip = (right panel) and 0.4 (left panel). On the z axis the variance 
er(p) of the Bose density p is plotted. For the MI phase the variance of the 2D density is zero, and the Bose charge has no 
fluctuation. 



The phase diagram in the p,, t plane is depicted in Fig[T] which is calculated for <p = and (p = 0.1. For the Mott 
insulator state the site occupation number p{nm) is equal to p and has integer values (see Mott lobes in Fig.[T] for 
p = 1, 2, 3). The variance of p is equal to zero for the MI state and has large fluctuations for the SF state (see the 
variance of p in Fig. [2]). 

The Bose-Hubbard model under a magnetic field has also been considered by Niemeyer, Freericks and Monien [24| . 
using strong coupling expansion. Strong coupling expansion utilizes a perturbative expansion in t/U, and is valid 
within the Mott insulating regime for small t/U and small values of the flux ip. The mean- field approach of this 
paper, on the other hand, is a self-consistent but uncontrolled approximation, which makes it possible to calculate 
physical quantities for both the insulating and the superfluid regimes. Although the two approximation methods 
have very different character, they generally yield qualitatively similar results. Very similar phase diagrams have been 
obtained by these two different approaches for the pure Bose-Hubbard model ([23| and [25|]), and the superlattice 
Bose-Hubbard model ( 26] and [27]). Indeed, we find that our mean-field treatment produces a qualitatively similar 
result to the strong coupling expansion for the Bose-Hubbard model under a magnetic field. However, our method 
allows us to extend the calculation to the full range of magnetic flux < ip < 1 , and calculate the superfluid density, 
and magnetization in the superfluid phase. 

We note that the transition point at zero magnetic flux is located at t c = 0.043 for p = 0.5 (c.f . Fig.[TJ) . This 
corresponds to U/At = 5.8 which is equal to the transition point determined by van Oosten et al. [281 ]. 

The magnetic flux breaks the temporal invariance of the system and destroys the coherence of the wave function in 
a 2D fermion system after a "flight" time proportional 29] to 1/B. For a Bose system, the magnetic flux can have the 
same effect over the SF coherent wave function, destroying the stability of the superfluid solution in the proximity of 
the transition point and leading to the SF-MI transition when ip is bigger than a critical value tp c . This is consistent 
with the increased area of the Mott lobes when the magnetic field is present as shown in Figs. [1] and O For all values 
of the magnetic flux we depicted the phase diagram in the i, ip plane for the different values of the chemical potential 
in Fig.d 

We note the above MI-SF transition for a constant t value, increasing the magnetic flux ip from zero to the critical 
transition point (p c . For small <p the optical lattice approximates a continuous system of charged bosons under the 
magnetic field and the superfluid-insulator transition is similar to the disappearance of the superconductivity when 
the external magnetic flux in a superconductor exceeding a critical value [3fJ. 

However, for the Bose-Hubbard Hamiltonian the discreteness of the lattice brings about other interesting features. 
The noninteracting spectrum [22| is periodic in <p with Atp — 1 and symmetric around the value <p — 0.5. We recover 
the same feature of the phase diagram in Fig.[3J 

One of the most important results of our calculation is displayed in Fig. [3] The oscillations of the critical hopping 
strength with changing magnetic field follow the oscillations of the bandwidth of the Hofstadter butterfly [2^ (See 
Fig.[4j. For instance, increasing ip the bandwidth of the Hofstadter butterfly system [13] is shrunk, meaning that the 
superfluid order appears at a smaller value of U (or higher value of t). For tp — > 0.5 the bandwidth of the Hofstadter 
butterfly is increased again and the critical value t c of the transition is increased. This suggests that in the phase 
diagram of Fig. [3] (calculated for U = 1), the critical point (U/t) c of the SF-MI transition is proportional to the 
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Phase Diagram at ^=0.3 and 0.4 
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FIG. 3: (color online) Phase diagram of Bose atoms confined in a 2D optical lattice in the t — ip plane for different values of the 
chemical potential, corresponding to the MI phase with p — 1 (top) and p = 2 (bottom). The curves correspond to the critical 
values t c of the SF-MI transition. For t > t c the system is in the superfluid phase and for t < t c in the MI phase. The phase 
diagram is periodic in the magnetic flux with Aip — 1 and symmetric around the flux value </? = 0.5. 

bandwidth of the Hofstadter butterfly [2^|. 

The loss of translational invariance of the charged bosons under a magnetic field is correlated with the appearance 
of the surface oscillations of the particle density in the superfluid phase. Although one expects a spontaneous breaking 
of the translational invariance in the superfluid phase, such as vortex states, in our mean-field approach this symmetry 
is explicitly broken by our gauge choice. Physically measurable quantities should be independent of the choice of the 
gauge, however our mean-field treatment yields spatially dependent parameters such as | ^ rl | 2 , which depend on the 
gauge choice. Still, it is not unreasonable to expect that spatially averaged quantities such as p s and magnetization 
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FIG. 4: Single particle energy bands as a function of magnetic flux per plaquette, the Hofstadter butterfly [22 



to be correctly captured by our approach. For certain values of the magnetic field we verified this expectation by 
using the symmetric gauge and corresponding square supercell. Also, as the insulator side of the transition is spatially 
uniform, explicit determination of the gauge should not strongly affect the Ml-superfluid phase boundaries. In Fig. [5] 
we show the variance of the on-site SF order parameter \I/ n when the hopping parameter t is increased, for different 
points of the lattice n. Even all the SF order parameters \l/ n exhibit the disappearance of the insulator order at 
the same critical value, t = t c , the rate of increase of ^> n for t > t c is different for different points of the lattice. 
There are points in the lattice that show a very low rate of increase of The figure clearly shows the surface 
oscillations of the 'charge' density in the optical lattice. It means that, for t > t c , there are regions in the lattice with 
small SF order parameter and small density fluctuations. It is possible to interpret these as the superfluid density 
oscillations due to the presence of vortices. However, close to the Mott transition where these regions are most 
pronounced, correlations between such regions may develop causing a phase transition which would not be captured 
by our mean-field treatment. 

The lower rate of increase of the superfluid order parameter vs. t, for t > t c , is also suggested in the scaled curves 
in Fig.Othat show the ratio p/p s versus t c /t. For t c /t ~< 1, the SF phase is strongly affected by the magnetic field 
presence, and we note the change of the curve slope in Fig. [S] for ip ^ 0. For t c /t — > 0, the magnetic field has no effect, 
and all Bose particles tend to condense (p/p s 1 for all ip values in Fig.[S]). We also notice that even a small magnetic 
field ip = 0.1 affects the superfluid density rather strongly. For such small magnetic fields our approach is less reliable, 
as the cell size used in our calculations is inversely proportional to the flux ip. While similar discontinuous effects at 
zero field have been discussed within the context of Josephson junction arrays [31] , we can not clarify the behavior 
near zero field due to the limitations of our numerical mean-field approach. 

Figure [7] shows the magnetization as a function of the hopping parameter t for different values of the magnetic 
flux (p = 0.1, 0.2, 0.3 and = 0.4 at the chemical potential p, — 0.5. We note that our approach is less reliable for 
the case of ip — 0.1, as the number of particles per site increase most quickly in this case. For t < i c , in the MI 
phase, the magnetization is equal to zero and the particle motion is frozen. For t > t c , in the SF phase, the curves 
exhibit negative values of magnetization, meaning the existence of persistent current flow. We note that the sign of 
the magnetization is related to the slope of the energy curve vs. magnetic flux in Hofstadter butterfly [3~3 |. 

The change of slope of the eigenstates as a function of <p in the Hofstadter butterfly gives rise to changes in the 
sign of the magnetization near the special values of the magnetic flux (for instance ip = 1/q; see Ref. [HI)- We do not 
expect this fine effect to be observable in the low ip limit. In this case, the small fluctuations of the <p for the real 
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FIG. 5: (color online) SF order parameter ^ n for different lattice sites at MI-SF transition. The chemical potential is p = 0.5 
corresponding to MI with p = 1 and /i = 1.5 corresponding to MI with p — 2. The blue line with crosses represents the 
surface average value For t > t c the second derivative d 2 ^/dt 2 changes sign for the points n = 4, 5. The magnetic flux is 
ip — 0.1 and the superfiuid order parameter has the periodicity >P„ = <&n+io. For the SF phase (at t > t c ) the lattice sites 
integer * 10 + 4 and integer * 10 + 5 exhibit very low values of the SF order parameter and low variance of p. 



system gives rise to the smeared graph of the spectrum [HJ . 

For <p S [0.45, 0.5] the magnetization can have an opposite sign (compared with the values for <p G [0, 0.45]) as the 
slope of the Hofstadter spectrum vs. ip clearly changes (see [22| or Fig. [3]). 
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FIG. 6: (color online) ps/p as a function of t c /t for different values of the magnetic flux. For ip = our curve is similar to 
Fig. 1 of Ref. [23|]. The change of slope at t ~> t c for p 7^ is related to the surface oscillations of the charged bosons. When 
the magnetic field is present, the superfluid phase exhibits surface region where the quantum fluctuations of SF density and 
the order parameter are small (compared to nearby regions; see Fig. [5} meaning that the local phase is closer to an insulator. 
It exhibit a lower ratio p 3 /p compared with the case ip = 0. 



H=0.5 

0.25 1 . . 




FIG. 7: (color online) Magnetization as a function of t for ip = 0.1, 0.2, 0.3, 0.4 (from top to bottom). In the MI phase, for 
t < t c , the magnetization is zero. Our numerical approach is more reliable at lower values of t, and not too small cp, where our 
computational supercell size is small and the number fluctuations at each lattice site are limited. MI and superfluid phases 
and t c are identified for ip = 0.1 on the figure. The value of t c for other ip can be identified as the point where magnetization 
becomes non-zero on the corresponding curve. 
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IV. SUMMARY 

We calculated the mean-field phase diagram of the 2D Bose-Hubbard Hamiltonian under a perpendicular magnetic 
field. The Mott insulator-superfluid transition is strongly affected, and the features around the transition point 
resemble the interesting property of the Hofstadter butterfly. The energy spectrum for noninteracting case exhibits 
periodicity with Atp = 1, symmetry around the value ip = integer/2 and striking oscillations [22j that lead to similar 
features of the MI-SF transition when the magnetic field is varied (see the phase diagram in t, ip plane in Fig. [3]). 

In the superfluid phase, at zero magnetic field the system has time translational invariance and the net local current 
of the Bose particles is zero because the reversed Bose paths are equally probable. Even for a small value of ip the time 
invariance is suddenly broken and the reversal paths of the coherent Bose atoms are not equally probable anymore (in 
fact only one of the two reversal paths is permitted; the other one corresponds to change of the sign of the magnetic 
flux). The persistent currents of the Bose particles appear, leading to nonzero value of the orbital magnetization. We 
note also the surface charge oscillations, that are not present in the SF phase with translational invariance when the 
magnetic field is set to zero. 

The surface charge oscillations are not present for MI phase. At low values of ip the SF phase exhibit negative 
susceptibility when the magnetic field is applied and at a critical value the system is driven into the insulating phase. 
This is the continuous limit of the model and the behavior of the Bose system shows features resembling to that 
of the Meissner effect in superconductors. For higher kinetic energy (large values of t/U) this scenario is not valid 
anymore, the magnetic flux loses its effect over the quantum state of the Bose system and the system preserves 
its SF characteristics, albeit with spatial oscillations of its order paremeter. We hope our work stimulates further 
experimental and theoretical interest in this model. 
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